home *** CD-ROM | disk | FTP | other *** search
/ IRIX 6.2 Development Libraries / SGI IRIX 6.2 Development Libraries.iso / dist / complib.idb / usr / share / catman / p_man / cat3 / complib / clatps.z / clatps
Text File  |  1996-03-14  |  7KB  |  199 lines
  1.  
  2.  
  3.  
  4. CCCCLLLLAAAATTTTPPPPSSSS((((3333FFFF))))                                                          CCCCLLLLAAAATTTTPPPPSSSS((((3333FFFF))))
  5.  
  6.  
  7.  
  8. NNNNAAAAMMMMEEEE
  9.      CLATPS - solve one of the triangular systems   A * x = s*b, A**T * x =
  10.      s*b, or A**H * x = s*b,
  11.  
  12. SSSSYYYYNNNNOOOOPPPPSSSSIIIISSSS
  13.      SUBROUTINE CLATPS( UPLO, TRANS, DIAG, NORMIN, N, AP, X, SCALE, CNORM,
  14.                         INFO )
  15.  
  16.          CHARACTER      DIAG, NORMIN, TRANS, UPLO
  17.  
  18.          INTEGER        INFO, N
  19.  
  20.          REAL           SCALE
  21.  
  22.          REAL           CNORM( * )
  23.  
  24.          COMPLEX        AP( * ), X( * )
  25.  
  26. PPPPUUUURRRRPPPPOOOOSSSSEEEE
  27.      CLATPS solves one of the triangular systems
  28.  
  29.      with scaling to prevent overflow, where A is an upper or lower triangular
  30.      matrix stored in packed form.  Here A**T denotes the transpose of A, A**H
  31.      denotes the conjugate transpose of A, x and b are n-element vectors, and
  32.      s is a scaling factor, usually less than or equal to 1, chosen so that
  33.      the components of x will be less than the overflow threshold.  If the
  34.      unscaled problem will not cause overflow, the Level 2 BLAS routine CTPSV
  35.      is called. If the matrix A is singular (A(j,j) = 0 for some j), then s is
  36.      set to 0 and a non-trivial solution to A*x = 0 is returned.
  37.  
  38.  
  39. AAAARRRRGGGGUUUUMMMMEEEENNNNTTTTSSSS
  40.      UPLO    (input) CHARACTER*1
  41.              Specifies whether the matrix A is upper or lower triangular.  =
  42.              'U':  Upper triangular
  43.              = 'L':  Lower triangular
  44.  
  45.      TRANS   (input) CHARACTER*1
  46.              Specifies the operation applied to A.  = 'N':  Solve A * x = s*b
  47.              (No transpose)
  48.              = 'T':  Solve A**T * x = s*b  (Transpose)
  49.              = 'C':  Solve A**H * x = s*b  (Conjugate transpose)
  50.  
  51.      DIAG    (input) CHARACTER*1
  52.              Specifies whether or not the matrix A is unit triangular.  = 'N':
  53.              Non-unit triangular
  54.              = 'U':  Unit triangular
  55.  
  56.      NORMIN  (input) CHARACTER*1
  57.              Specifies whether CNORM has been set or not.  = 'Y':  CNORM
  58.              contains the column norms on entry
  59.              = 'N':  CNORM is not set on entry.  On exit, the norms will be
  60.  
  61.  
  62.  
  63.                                                                         PPPPaaaaggggeeee 1111
  64.  
  65.  
  66.  
  67.  
  68.  
  69.  
  70. CCCCLLLLAAAATTTTPPPPSSSS((((3333FFFF))))                                                          CCCCLLLLAAAATTTTPPPPSSSS((((3333FFFF))))
  71.  
  72.  
  73.  
  74.              computed and stored in CNORM.
  75.  
  76.      N       (input) INTEGER
  77.              The order of the matrix A.  N >= 0.
  78.  
  79.      AP      (input) COMPLEX array, dimension (N*(N+1)/2)
  80.              The upper or lower triangular matrix A, packed columnwise in a
  81.              linear array.  The j-th column of A is stored in the array AP as
  82.              follows:  if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
  83.              if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
  84.  
  85.      X       (input/output) COMPLEX array, dimension (N)
  86.              On entry, the right hand side b of the triangular system.  On
  87.              exit, X is overwritten by the solution vector x.
  88.  
  89.      SCALE   (output) REAL
  90.              The scaling factor s for the triangular system A * x = s*b,  A**T
  91.              * x = s*b,  or  A**H * x = s*b.  If SCALE = 0, the matrix A is
  92.              singular or badly scaled, and the vector x is an exact or
  93.              approximate solution to A*x = 0.
  94.  
  95.      CNORM   (input or output) REAL array, dimension (N)
  96.  
  97.              If NORMIN = 'Y', CNORM is an input argument and CNORM(j) contains
  98.              the norm of the off-diagonal part of the j-th column of A.  If
  99.              TRANS = 'N', CNORM(j) must be greater than or equal to the
  100.              infinity-norm, and if TRANS = 'T' or 'C', CNORM(j) must be
  101.              greater than or equal to the 1-norm.
  102.  
  103.              If NORMIN = 'N', CNORM is an output argument and CNORM(j) returns
  104.              the 1-norm of the offdiagonal part of the j-th column of A.
  105.  
  106.      INFO    (output) INTEGER
  107.              = 0:  successful exit
  108.              < 0:  if INFO = -k, the k-th argument had an illegal value
  109.  
  110. FFFFUUUURRRRTTTTHHHHEEEERRRR DDDDEEEETTTTAAAAIIIILLLLSSSS
  111.      A rough bound on x is computed; if that is less than overflow, CTPSV is
  112.      called, otherwise, specific code is used which checks for possible
  113.      overflow or divide-by-zero at every operation.
  114.  
  115.      A columnwise scheme is used for solving A*x = b.  The basic algorithm if
  116.      A is lower triangular is
  117.  
  118.           x[1:n] := b[1:n]
  119.           for j = 1, ..., n
  120.                x(j) := x(j) / A(j,j)
  121.                x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j]
  122.           end
  123.  
  124.      Define bounds on the components of x after j iterations of the loop:
  125.         M(j) = bound on x[1:j]
  126.  
  127.  
  128.  
  129.                                                                         PPPPaaaaggggeeee 2222
  130.  
  131.  
  132.  
  133.  
  134.  
  135.  
  136. CCCCLLLLAAAATTTTPPPPSSSS((((3333FFFF))))                                                          CCCCLLLLAAAATTTTPPPPSSSS((((3333FFFF))))
  137.  
  138.  
  139.  
  140.         G(j) = bound on x[j+1:n]
  141.      Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}.
  142.  
  143.      Then for iteration j+1 we have
  144.         M(j+1) <= G(j) / | A(j+1,j+1) |
  145.         G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] |
  146.                <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | )
  147.  
  148.      where CNORM(j+1) is greater than or equal to the infinity-norm of column
  149.      j+1 of A, not counting the diagonal.  Hence
  150.  
  151.         G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | )
  152.                      1<=i<=j
  153.      and
  154.  
  155.         |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| )
  156.                                       1<=i< j
  157.  
  158.      Since |x(j)| <= M(j), we use the Level 2 BLAS routine CTPSV if the
  159.      reciprocal of the largest M(j), j=1,..,n, is larger than
  160.      max(underflow, 1/overflow).
  161.  
  162.      The bound on x(j) is also used to determine when a step in the columnwise
  163.      method can be performed without fear of overflow.  If the computed bound
  164.      is greater than a large constant, x is scaled to prevent overflow, but if
  165.      the bound overflows, x is set to 0, x(j) to 1, and scale to 0, and a
  166.      non-trivial solution to A*x = 0 is found.
  167.  
  168.      Similarly, a row-wise scheme is used to solve A**T *x = b  or A**H *x =
  169.      b.  The basic algorithm for A upper triangular is
  170.  
  171.           for j = 1, ..., n
  172.                x(j) := ( b(j) - A[1:j-1,j]' * x[1:j-1] ) / A(j,j)
  173.           end
  174.  
  175.      We simultaneously compute two bounds
  176.           G(j) = bound on ( b(i) - A[1:i-1,i]' * x[1:i-1] ), 1<=i<=j
  177.           M(j) = bound on x(i), 1<=i<=j
  178.  
  179.      The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we add
  180.      the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1.  Then the
  181.      bound on x(j) is
  182.  
  183.           M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) |
  184.  
  185.                <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| )
  186.                          1<=i<=j
  187.  
  188.      and we can safely call CTPSV if 1/M(n) and 1/G(n) are both greater than
  189.      max(underflow, 1/overflow).
  190.  
  191.  
  192.  
  193.  
  194.  
  195.                                                                         PPPPaaaaggggeeee 3333
  196.  
  197.  
  198.  
  199.